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We describe recent results concerning the behavior of lattice QCD with light dynamical Wilson and Staggered 
quarks. We show that it is possible to reach regions of parameter space with light pions PS 0.2/a using Wilson 
fermions. If the Hybrid Molecular Dynamics (HMD) algorithm is used with the same parameters it gives incorrect 
results. We also present preliminary results using a higher-order integration scheme. 



1. Introduction 

In this paper, we set out with the objective of 
examining the effects of 6H errors in the integra- 
tion of Hamilton's equations of motion. We have 
performed tests using 2 flavors of Wilson fermions 
at /3 = 5.30, k = 0.1675 and 0.1677 on a 16 3 x 32 
lattice. We have also studied four flavors of Stag- 
gered fermions at [3 = 5.05, m q = 0.005 on an 
8 3 x 16 lattice. Previous studies of step-size er- 
rors including fermions were limited to relatively 
small lattices and large pion masses. Our goal is 
to extend these studies to parameter regions that 
will be realistic on future machines. 

The standard methods of simulating QCD 
with dynamical fermions are Hybrid Monte Carlo 
(HMC) via the (^-algorithm [1] and Hybrid Molec- 
ular Dynamics (HMD) via the i?-algorithm [2]. 
The major distinction between the two algo- 
rithms is how the fermionic determinant is simu- 
lated and how the configuration obtained by ap- 
proximately integrating Hamilton's equations of 
motion is handled. The nomenclature is some- 
what confusing. The expression Hybrid Molec- 
ular Dynamics is principally used to name the 
i?-algorithm; however, it is also used to describe 
the (^-algorithm without the acceptance test. The 
principle advantage of the (^-algorithm with the 
acceptance test (HMC) is that the algorithm is 
exact - no finite step-size errors are introduced 
into the equilibrium probability distribution. 

Our tests are limited to the Hybrid Monte 
Carlo algorithm where we can compare the al- 



* Speaker at the conference 



gorithm with and without the Metropolis accep- 
tance. We will call the algorithm without the 
Metropolis acceptance the HMD algorithm. All 
our computations were done on the Connection 
Machine CM-2 at SCRI using a conjugate gra- 
dient linear equation solver (written in CMIS) 
which runs at > 5 Gigaflops. 

A major difficulty in simulating QCD is the 
long exponential autocorrelation times for the 
system to decay into equilibrium and the long 
integrated autocorrelation times while in equilib- 
rium. In fact, our results for Wilson fermions with 
the lightest pion [m w ps 0.2/a) indicate a decay 
time of about 700 unit-length trajectories. We 
should expect this behavior. The pion is the light- 
est excitation of QCD, because it is trying to be 
a Goldstone boson. Unlike the quenched approx- 
imation where the lightest dynamical excitation 
is a glueball, for full QCD the pion mass is the 
inverse correlation length, m w = l/£. Thus, for 
= 0.2/a we expect that the correlation length 
£ ps 5a in equilibrium; it is only reasonable to ex- 
pect an algorithm to take much longer to develop 
such long distance correlations than it takes to 
equilibrate a quenched system with £ ps a. The- 
oretically, we expect the characteristic relaxation 
time to be the same as the autocorrelation time, 
which in turn is proportional to £ z where z is the 
dynamical critical exponent. For HMC it is be- 
lieved that 1 < z < 2 [5]. 

How can we measure the correlation length 
of an ensemble of configurations? The obvious 
method of measuring £ is to measure the pion 
mass by measuring the asymptotic decay of a pion 
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correlation function (ir(0)ir(t)) . If we have a set 
of equilibrated configurations this is a completely 
valid procedure, but consider what would happen 
if the Markov process had not yet converged to 
the true equilibrium distribution. Because the 
system starts from a configuration which does 
not have the true long distance correlations built 
in, the unequilibrated configurations will tend to 
have £ < l/m w . Roughly speaking we can say 
that they correspond more closely to the equilib- 
rium distribution of QCD with a larger value of 
the quark mass m q — we shall call this the "sea" 
quark mass - than the actual value appearing in 
the action. Of course, the actual distribution of 
some set of unequilibrated configurations is some 
complicated mess, but it is certainly reasonable 
to assume that such a shift in m q is the dominant 
effect when near m cr8t = for Staggered fermions 
or k c in the case of Wilson fermions. 

The quark mass does not only appear in the 
action, it is also explicitly present in the pion op- 
erator 7r used to measure the pion mass. Since 
this quark mass has nothing to do with the dy- 
namics of the system we shall call it the "valence" 
quark mass. For our computations we expect that 
the sea m q approaches the valence m q from above 
as the system approaches equilibrium: in other 
words the correlation length estimated using the 
correlation function for a "valence" pion will be 
larger than the actual correlation length of the 
system until the system reaches equilibrium. 

For a quenched computation or for an inexact 
dynamical quark algorithm valence and sea pa- 
rameters are different even in equilibrium. One 
illustration of these ideas is that it is easy to mea- 
sure a light pion mass on a quenched configura- 
tion, even though the correlation length is much 
smaller than the inverse pion mass. 

The existence of a light valence pion is not ev- 
idence for a system containing light dynamical 
pions unless the system can be shown to have 
equilibrated to the correct probability distribu- 
tion. A good indicator we have found is blocked 
gauge field operators (using Teper blocking), for 
example the blocked spatial plaquette. Blocked 
plaquettes do not need an extra parameter and 
are sensitive to the longer range properties of the 
system. The effects of dynamical quarks should 



appear as a modification of the quenched QCD 
beta function at longer distances. The blocked 
observables indicate coherence and correlation in 
the system. 

The proceeding general arguments are applica- 
ble to both Wilson and Staggered fermions. The 
difference between the two fermions is at which 
quark mass scale should the step-size errors sig- 
nificantly effect the equilibrium probability distri- 
bution. The fermionic contribution to the molec- 
ular dynamics force is 

X (Aft M)" 1 — [Aft Af] (Aft Af)" 1 X (1) 

where Af = M(U(t)), U(t) is the gauge field 
configuration at integration time t and x is the 
pseudo-fermion field. When the quark mass ap- 
proaches criticality the fermionic matrix becomes 
singular and the smallest eigenvalue dominates 
the force. The lowest eigenvalue Ao of Aft Af is 
related to the pion mass, hence we argue that the 
step-size errors in the integration are amplified by 
some function of l/m w . Wilson fermions do not 
have chiral symmetry hence Ao can wander arbi- 
trarily close to zero during integration. On the 
other hand, for Staggered fermions Ao > for all 
configurations with m q > 0. Most likely the func- 
tional dependence on l/m w is different for Wilson 
and Staggered fermions, but in any case we ex- 
pect very small pion masses and a small sea quark 
mass will be required for Staggered fermions be- 
fore the magnified step-size errors will be an issue. 

For our Wilson fermion tests, the first exercise 
we had to undertake was to find a suitable value 
for k. We ran at k = 0.1680, 0.1677, 0.1675, and 
0.1670: our results indicate that m w < 0.2/a for 
the first two values, indicating that k c is prob- 
ably in this region. If we extrapolate recently 
published results [3] we find that k c ps 0.1676, in 
rough agreement with our data. For k = 0.1670, 
0.1675 and 0.1677 we find that « 0.47/a, 
0.33/a and 0.2/a, respectively. 

In Figure 1 we show the pion mass and ac- 
ceptance rate as a function of MD time for our 
run at k = 0.1677. Here, Rmd is the accuracy 
of the normalized residual vector during integra- 
tion. It is clear the system has a large exponen- 
tial autocorrelation time and probably has only 
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Figure 1. Time history of the effective pion mass 
(distance 8) for [3 = 5.30, k = 0.1677 and St = 
0.0069 on a 16 3 x 32 lattice. Symbols are: (x) 
HMC at Rmd = 10 -5 before trajectory 606 and 
R MD = 3 x 10" 7 after 606, (+) HMD at R M d = 
10" 5 and (□) HMD at R M d = 3 x 10" 7 . Also 
shown is the blocksize= 40 (O) HMC acceptance 
uncorrected for autocorrelations. 



just reached equilibrium. It is also clear that as 
the correlation length grows the HMC acceptance 
rate falls, which requires reducing the integration 
step-size and/or the conjugate gradient accuracy 
in order to keep a reasonable acceptance rate. A 
preliminary study of the system indicates that 
this is primarily due to the increasing effect of 
high frequency fluctuations in the fermionic con- 
tributions to the action. In all our Wilson runs, 
we used a linear extrapolation of the proceeding 
two solutions for the initial guess of the conjugate 
gradient inverter. 

We decided it would be useful to see what re- 
sults we would get if we used the HMD algo- 
rithm instead of HMC for this system (i.e., if 
we omitted the Monte Carlo accept/reject step 
for Rmd = 3 x 10 -7 ). It is commonly accepted 
that this should introduce small errors in phys- 
ical quantities of order St 2 ; to our surprise the 
HMD algorithm with the same parameters gave 



wrong answers for the pion mass. We also per- 
formed a HMD test using Rmd = 10 -5 and found 
the algorithm gave completely wrong answers for 
the pion mass. In fact other hints that the sys- 
tem contained light pions, such as the large num- 
ber of CG iterations required per step to reach 
a given residual, also rapidly disappeared when 
running the HMD algorithm at R M d = 10" 5 . No 
shown is a test of when we reduced the step size 
to St = 0.004 at R M d = 10" 5 , the HMD results 
were virtually identical to St = 0.0069 indicating 
the SH errors were dominated by the CG accu- 
racy. We haven't performed any St HMD tests at 
Rmd = 3 x 10 -7 ; however, we expect the CG ac- 
curacy to muddy the results. An even clearer in- 
dicator of SH errors introduced is via the blocked 
spatial plaquette. At k = 0.1675 and 0.1677 and 
using Rmd = 3 x 10 -7 , the blocked plaquette for 
HMD deviated systematically from HMC when 
the acceptance test was turned off while the pion 
mass for k = 0.1675 with HMD deviated only 
slightly from HMC. 

These tests show that St, CG accuracy and the 
acceptance are all very crucial for the algorithm 
to give correct physics. Tests of only step-size 
extrapolations are not complete enough since the 
CG accuracy is important for not just the mini- 
mization of SH errors but for reversibility as well. 

The parameters used in our Staggered simula- 
tions were based on the work of the MTc collabo- 
ration [4]. We used an 8 3 x 16 lattice at (3 = 5.05 
and m q = 0.005 and observed m w = 0.188(3) and 
m p k, 1.0. For HMC, we used R M d = 5 x 10" 7 
and a linear extrapolation of the initial guess. To 
emphasize the SH errors we used Rmd = 5x 10 -4 
for HMD with a zero initial guess to maintain 
reversibility. We didn't find any statistically sig- 
nificant effects between HMD and HMC after 500 
and 1000 trajectories each, respectively. The pion 
mass here is too small and the rho is too large 
to be realistic of the continuum limit, however. 
With this small lattice, we probably shouldn't ex- 
pect the sea quark mass to be very small. We 
are currently running at the MTc parameters of 
j3 = 5.35 and m q < 0.01 on a 16 3 x 32 lattice 
which we hope will provide a test more useful for 
calculations on future machines. 

In a separate study, we carried out a more accu- 
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rate MD integration with Wilson fermions using 
the higher-order integration schemes [6], with the 
hope of finding a cheap way of increasing the ac- 
ceptance rate for HMC. Theoretical analyses of 
these algorithms were based upon the idea that 
because the change in energy over a trajectory 
would be 6H = 0(8t 4 ) instead of 6H = 0(8t 2 ) 
one could use a much larger step size to get the 
same acceptance rate. This does not work in 
practice: while the Campostrini scheme does in- 
deed give much smaller values for 8H, the inte- 
gration becomes unstable if any attempt is made 
to increase St. Indeed, the Campostrini method 
proved unstable unless the step size was reduced 
slightly from that used in our HMC runs. This 
is, of course, the behavior we should have ex- 
pected. The limit on the step size is given by 
the constraint from stability (u>6t <C 1 where u> 
is the highest frequency of the physical system) 
and not by the erfc(const ■ 8t' 2 \/V) dependence 
of the acceptance rate on the volume (at least for 
the volumes we have used). 

Initially we thought the extra computation re- 
quired for the Campostrini method was small be- 
cause the interpolated solutions of the previous 
CG solutions were very accurate. This was erro- 
neous, as we had overlooked the fact that the del- 
icate cancellations required by the Campostrini 
method required a much more accurate measure- 
ment of H for each step, and thus a much lower 
CG residual. When we lowered the CG residual to 
the value required to make the CG solution error 
contribution to H smaller than the true step-size 
errors, we found that the magnitude of 6H was 
greatly reduced but at the cost of many more CG 
iterations. Sadly we found that unless enough ex- 
tra CG iterations were performed the use of inter- 
polated initial guesses for the CG solver induced 
an irreversibility in the HMC algorithm which led 
to clearly incorrect results. 

Our conclusions on this issue are that one must 
either make sure that the CG residual is small 
enough to avoid such systematic biases or one 
must use a time-symmetric initial guess for the 
CG solution such as zero. Both options make the 
Campostrini method much more expensive than 
simple leapfrog integration. Which method is ul- 
timately cheaper is still under investigation. 



If we had carried out an HMD calculation with 
similar 6H errors used in the present HMC com- 
putations, we would have been forced to use a 
larger value of k in order to measure a light 
pion mass. This mass would have been much 
smaller than the inverse correlation length actu- 
ally present in the configurations. Indeed, we 
may well have been measuring valence observ- 
ables on an almost-quenched system, so we would 
have concluded that physics with light dynamical 
quarks was very similar to quenched physics with 
the same masses. We would also have found a 
much smaller relaxation time and autocorrelation 
time, and thus been misled into thinking that full 
QCD calculation were much cheaper than they 
really are. 

What are the implications of this for staggered 
quarks calculations done using the HMD algo- 
rithm? While our Wilson quark results do not 
tell us directly about the behavior of staggered 
quarks — the effects we see might possibly just 
be artifacts of Wilson quark dynamics — they 
lead us to suggest that it is necessary for a care- 
ful zero step-size extrapolation to be done for any 
HMD calculation, with special attention required 
to verify that the system is truly in equilibrium. 
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